###############
# Replication Code for: The Electoral Consequences of International Migration in Sending Countries
###############

# Load libraries 
if (!require("pacman")) install.packages("pacman", repos = "http://cran.us.r-project.org")
pacman::p_load(stargazer,tidyverse,tidyr,plm, sandwich, lmtest, ggpubr)

#Working directory
argv <- getwd()
setwd(argv)

###Data
#Figure 1, Table 1-2
lits<- read.csv("lits_2_selected.csv", stringsAsFactors = F)
ess<- read.csv("ess_5_9_selected.csv", stringsAsFactors = F)

#Table 4-5
ches_emig_dat<- read.csv("ches_emig_dat.csv", stringsAsFactors = F)

#Table 5
pl_sample <- read.csv("pl_iv.csv", stringsAsFactors = F)

#Table 6
polpan<- read.csv("polpan_clean.csv", stringsAsFactors = F)

###########################################
#Figure 1
###########################################
#Figure 1 (LiTs)
lits<- lits %>% filter(!is.na(emig)) %>%
  mutate(emig2=ifelse(emig==1, "Willing", "Not Willing"))


#common theme for the figures
mytheme = list(
  theme_classic()+
    theme(legend.position =  "none",
          axis.text.x = element_text(size=27),
          axis.text.y = element_text(size=27, angle=45), 
          axis.title.x = element_blank(), axis.title.y = element_blank()))


#Age Density Plot
fig1a1<- lits %>% filter(!is.na(age)) %>% ggplot()+
  geom_density(aes(x=as.numeric(age), group = as.factor(emig2), color = as.factor(emig2),
    fill=as.factor(emig2)), alpha=0.4, adjust=1, na.rm = T)+
  scale_fill_manual(values=c("gray","red"))+
  scale_colour_manual(values=c("darkgray", "lightpink"))+mytheme+
  scale_x_continuous(breaks=seq(18,99,20))+
  scale_y_continuous(breaks=seq(0, 0.06, 0.01))


#Education Density Plot
educ.label =c("Primary", "", "", "","",  "BA", "MA+")

fig1b1<- lits %>% filter(!is.na(educ)) %>% 
  ggplot()+ geom_density(aes(x=factor(educ), 
            group = as.factor(emig2), color = as.factor(emig2), fill=as.factor(emig2)),
      adjust=3,alpha=0.4,na.rm = T)+mytheme+
  scale_fill_manual(values=c("darkgray","red"))+
  scale_colour_manual(values=c("darkgray","lightpink"))+ylim(c(0, 0.4))+
scale_x_discrete(labels= educ.label)

anti_imm_label <- c("Pro-", "Neutral", "Anti-")

fig1c1 <- lits %>% filter(!is.na(anti_imm)) %>%
  ggplot()+geom_density(aes(x=factor(anti_imm), 
  group = as.factor(emig2),color = as.factor(emig2), fill=as.factor(emig2)), 
  adjust=4,alpha=0.4,na.rm = T)+
  scale_fill_manual(values=c("gray", "red"))+
  scale_colour_manual(values=c("darkgray", "lightpink"))+mytheme+
  scale_x_discrete(labels=anti_imm_label)+ylim(0, 0.5)

########Figure 1 (ESS)

#Age Density Plot
fig1a2<- ess %>% filter(!is.na(forborn) & !is.na(agem_correct)) %>%
  mutate(forborn=ifelse(forborn==0,  "Natives","Emigrants")) %>%
  mutate(forborn=factor(forborn, levels = c("Natives", "Emigrants"))) %>%
  ggplot()+geom_density(aes(x=agem_correct, group = forborn, color =forborn, fill=forborn), 
  adjust=3, alpha=0.4,na.rm = T)+mytheme+
  scale_fill_manual(values=c("gray", "red"))+
  scale_colour_manual(values=c("darkgray","lightpink"))+
  scale_x_continuous(breaks=seq(18,99,20))+scale_y_continuous(breaks = seq(0, 0.04, 0.01))

#Education Density Plot
fig1b2<- ess %>% filter(!is.na(forborn) & !is.na(eisced)) %>% filter(eisced !=0) %>% 
  mutate(forborn=ifelse(forborn==1, "Emigrants", "Natives")) %>%
  mutate(forborn=factor(forborn, levels = c("Natives", "Emigrants"))) %>%
  ggplot()+geom_density(aes(x=as.factor(eisced), group = forborn, color = forborn, fill=forborn),
  adjust=3.5,alpha=0.4,na.rm = T)+
  scale_fill_manual(values=c("gray", "red"))+
  scale_colour_manual(values=c("darkgray","lightpink"))+
  scale_y_continuous(breaks=seq(0, 0.4, 0.1))+
  scale_x_discrete(label=c("Primary", "", "", "", "", "BA", "MA+"))+mytheme

#Anti-Immigration Density Plot
anti_imm_label = c("Pro", "", "", "Anti")

fig1c2<- ess %>% filter(!is.na(forborn) & !is.na(anti_race_imm)) %>%
  mutate(forborn=ifelse(forborn==1, "Emigrants", "Stayers")) %>%
  mutate(forborn=factor(forborn, levels = c("Stayers", "Emigrants"))) %>%
  ggplot()+geom_density(aes(x=factor(anti_race_imm), 
        group = forborn, color = forborn, fill=forborn), adjust=5,alpha=0.4,na.rm = T)+
  scale_fill_manual(values=c("gray", "red"))+
  scale_colour_manual(values=c("darkgray","lightpink"))+mytheme+
  scale_x_discrete(labels=anti_imm_label)

###Table 1
#Table 1, model 1
t1_m1<- glm(emig ~ age+Female+educ+unemployed+sat_eco+
              as.factor(countryname), data = lits, family = "binomial")

#Table 1, model 2
t1_m2<- glm(emig ~ age+Female+educ+unemployed+
              sat_eco+anti_imm+support_dem+rlg_mem+vote+
              as.factor(countryname), data = lits, family = "binomial")

#Table 1: The Characteristics of (Potential) Emigrants (LiTs)
t1<- stargazer::stargazer(t1_m1, t1_m2, 
    omit = c('countryname', 'Constant'), omit.stat = 'aic', 
    add.lines = list(c("Country FE", "\\checkmark", "\\checkmark")), 
    covariate.labels = c("Age", "Female", "Education", "Unemployed", "Satisfied w Econ", 
        "Anti-Immigrant", "Support for Democracy", "Religiosity","Vote"), 
    dep.var.caption = "Willing to Emigrate", dep.var.labels.include = F)

########################################################
#Table 2: The Characteristics of Emigrants (ESS Wave 5-9)
t2_glm1<- glm(formula = forborn ~ agem_correct+female+eisced+anti_race_imm+lrscale+rlgdgr+as.factor(cntry2), data = ess, family = "binomial")

t2_glm2<- glm(formula = forborn ~ agem_correct+female+anti_race_imm+eisced+lrscale+rlgdgr+as.factor(Year)+as.factor(cntry2), data = ess, family = "binomial")

t2<- stargazer::stargazer(t2_glm1, t2_glm2,  omit = c("cntry2", "Constant", "Year"), 
omit.stat = c("f", "ser", "aic"), covariate.labels = c("Age (of arrival)", "Female",
        "Education", "Anti-Immigrant", "Ideology", "Religiosity"),
dep.var.caption = "Dependent variable: Emigrants", dep.var.labels.include = F,
add.lines = list(c("Country (of origin) FE", "\\checkmark", "\\checkmark"), c("Year FE", "", "\\checkmark")))

########Table 4

t4_fr0 <- lm(data = ches_emig_dat %>% filter(CHES == "rad right"), voteshare ~ emig.lag.1.voters+as.factor(nuts_code2013)+as.factor(Year))
t4_se0<- coeftest(t4_fr0, vcov=vcovHC(t4_fr0, type="HC1"))

t4_fr1 <- lm(data = ches_emig_dat %>% filter(CHES == "rad right"), voteshare ~ emig.lag.1.voters+imig.lag.1.voters+gdp.lag.1+log(Pop)+unemp.percent.lag.1.n2+c_tf_r.lag.1+as.factor(nuts_code2013)+as.factor(Year))
t4_se1<- coeftest(t4_fr1, vcov=vcovHC(t4_fr1, type="HC1"))

t4_fr2 <- lm(data = ches_emig_dat %>% filter(CHES == "rad right"), voteshare ~ emig.lag.1.voters+imig.lag.1.voters+gdp.lag.1+log(Pop)+unemp.percent.lag.1.n2+c_tf_r.lag.1+voteshare.w1+as.factor(nuts_code2013)+as.factor(Year))
t4_se2<- coeftest(t4_fr2, vcov=vcovHC(t4_fr2, type="HC1"))

t4_fr3 <- lm(data = ches_emig_dat %>% filter(CHES == "rad right"), voteshare ~ emig.lag.1.voters+imig.lag.1.voters+gdp.lag.1+log(Pop)+unemp.percent.lag.1.n2+c_tf_r.lag.1+voteshare.w1+WB_rem.lag.1+as.factor(nuts_code2013)+as.factor(Year))
t4_se3<- coeftest(t4_fr3, vcov=vcovHC(t4_fr3, type="HC1"))

#Table 4
t4<- stargazer::stargazer(t4_fr0, t4_fr1,  t4_fr2, t4_fr3, 
  omit=c("nuts_code2013","Year","Constant", "voteshare.w1"),omit.stat=c("f", "ser","adj.rsq", "rsq"), se = list(t4_se0[,2], t4_se1[,2], t4_se2[,2], t4_se3[,2]), align=TRUE,
  order = c("emig.lag.1", "imig.lag.1", "gdp.lag.1", "unemp.percent.lag.1.n2","Pop", "c_tf_r.lag.1","WB_rem.lag.1"), covariate.labels = c("Emigration", "Immigration","GDP","Unemployment", "Population","Current Transfers","(National) Remittances"), add.lines = list(c("NUTS FE", "\\checkmark", "\\checkmark", "\\checkmark", "\\checkmark"), c("Year FE", "\\checkmark", "\\checkmark", "\\checkmark", "\\checkmark"), c("Lagged DV", "", "","\\checkmark" ,"\\checkmark")))

######Table 5 

pl_sample<- left_join(ches_emig_dat %>% filter(cntry == "PL"), pl_sample)

t5_lm0<- lm(data=pl_sample %>% filter(CHES == "rad right") , voteshare ~ emig.lag.1.voters+as.factor(nuts_code2013)+as.factor(Year))
t5_lm0_se<- coeftest(t5_lm0, vcov=vcovHC(t5_lm0, type="HC1"))

t5_iv0 <- AER::ivreg(data=pl_sample %>% filter(CHES == "rad right"), voteshare ~ emig.lag.1.voters+as.factor(nuts_code2013)+as.factor(Year) |. - emig.lag.1.voters+as.factor(nuts_code2013)+as.factor(Year)+(unemp*emig.share.t1))
t5_iv0_diagnose<- summary(t5_iv0, diagnostics = T)
t5_iv0_se<- coeftest(t5_iv0, vcov=vcovHC(t5_iv0, type="HC1"))

t5_lm1<- lm(data=pl_sample %>% filter(CHES == "rad right") , voteshare ~ emig.lag.1.voters+imig.lag.1.voters+gdp.lag.1+unemp.percent.lag.1.n2+c_tf_r.lag.1+log(Pop)+as.factor(nuts_code2013)+as.factor(Year))
t5_lm1_se<- coeftest(t5_lm1, vcov=vcovHC(t5_lm1, type="HC1"))


t5_iv1 <- AER::ivreg(data=pl_sample %>% filter(CHES == "rad right"), voteshare ~ emig.lag.1.voters+imig.lag.1.voters+gdp.lag.1+unemp.percent.lag.1.n2+c_tf_r.lag.1+log(Pop)+as.factor(nuts_code2013)+as.factor(Year) -1 |. - emig.lag.1.voters+imig.lag.1.voters+gdp.lag.1+unemp.percent.lag.1.n2+c_tf_r.lag.1+log(Pop)+as.factor(nuts_code2013)+as.factor(Year)+(unemp*emig.share.t1))
t5_iv1_diagnose<- summary(t5_iv1, diagnostics = T)
t5_iv1_se<- coeftest(t5_iv1, vcov=vcovHC(t5_iv1, type="HC1"))


t5_lm2<- lm(data=pl_sample %>% filter(CHES == "rad right"), voteshare ~ emig.lag.1.voters+imig.lag.1.voters+gdp.lag.1+unemp.percent.lag.1.n2+voteshare.w1+c_tf_r.lag.1+log(Pop)+as.factor(nuts_code2013)+as.factor(Year))
t5_lm2_se<- coeftest(t5_lm2, vcov=vcovHC(t5_lm2, type="HC1"))


t5_iv2 <- AER::ivreg(data=pl_sample %>% filter(CHES == "rad right"), voteshare ~ emig.lag.1.voters+imig.lag.1.voters+gdp.lag.1+unemp.percent.lag.1.n2+c_tf_r.lag.1+log(Pop)+voteshare.w1+as.factor(nuts_code2013)+as.factor(Year) -1 |. - emig.lag.1.voters+imig.lag.1.voters+c_tf_r.lag.1+as.factor(nuts_code2013)+as.factor(Year)+gdp.lag.1+unemp.percent.lag.1.n2+voteshare.w1+(unemp*emig.share.t1))
t5_iv2_diagnose<- summary(t5_iv2, diagnostics = T)
t5_iv2_se<- coeftest(t5_iv2, vcov=vcovHC(t5_iv2, type="HC1"))

t5 <- stargazer::stargazer(t5_lm0, t5_iv0_se, t5_lm1, t5_iv1_se, t5_lm2, t5_iv2_se, 
  omit=c('nuts_code2013', 'Year', 'Constant', 'voteshare.w1'), 
  omit.stat = c("ser", "f", "rsq", "adj.rsq"), 
  se = list(t5_lm0_se[,2], t5_iv0_se[,2],t5_lm1_se[,2], t5_iv1_se[,2], t5_lm2_se[,2], t5_iv2_se[,2]), 
  dep.var.labels.include = F, model.names = F,add.lines = list(c("NUTS 3 FE", "\\checkmark", "\\checkmark", "\\checkmark","\\checkmark","\\checkmark", "\\checkmark"), c("Year FE", "\\checkmark", "\\checkmark", "\\checkmark","\\checkmark","\\checkmark", "\\checkmark"), c("Lagged DV", "", "", "", "","\\checkmark", "\\checkmark"), c("Estimator", "OLS", "IV", "OLS", "IV", "OLS", "IV"), c("First Stage F", "", round(t5_iv0_diagnose$diagnostics[1,3], 3),"", round(t5_iv1_diagnose$diagnostics[1,3], 3),"", round(t5_iv2_diagnose$diagnostics[1,3], 3))),covariate.labels = c("Emigration", "Immigration", "GDP", "Unemployment", "Current Transfers", "Population"))


########################################################
#Table 6

#Table 6, model 1
t6_within1 <- plm(data= polpan, model = "within", index=c("ANONID", "Year"),  curvote.pis1 ~ AGE+as.factor(GENDER)+BA_plus+nojob+emig.lag.1.percent, 
                  cluster='nuts_code2013')
t6_w1se<- coeftest(t6_within1, vcov=vcovHC(t6_within1, type="HC1"))

#Table 6, model 2
t6_within2 <- plm(data= polpan, model = "within", index=c("ANONID", "Year"), curvote.pis1 ~ AGE+as.factor(GENDER)+BA_plus+nojob+emig.lag.1.percent+imig.lag.1.percent, cluster='nuts_code2013')
t6_w2se<- coeftest(t6_within2, vcov=vcovHC(t6_within2, type="HC1")) 

#Table 6, model 3
t6_within3 <- plm(data = polpan, model = "within", index=c("ANONID", "Year"), 
 curvote.pis1 ~ AGE+BA_plus+nojob+emig.lag.1.percent+imig.lag.1.percent+gdp.lag.1.1000, cluster='nuts_code2013')
t6_w3se<- coeftest(t6_within3, vcov=vcovHC(t6_within3, type="HC1")) 


t6_within4 <- plm(data= polpan, model = "within", index=c("ANONID", "Year"), 
curvote.pis1 ~ AGE+BA_plus+nojob+emig.lag.1.percent+imig.lag.1.percent+gdp.lag.1.1000+unemp.percent.lag.1, cluster = 'nuts_code2013')
t6_w4se<- coeftest(t6_within4, vcov=vcovHC(t6_within4, type="HC1")) 


t6_within5 <- plm(data= polpan, model = "within", index=c("ANONID", "Year"), 
  curvote.pis1 ~ AGE+BA_plus+nojob+emig.lag.1.percent+imig.lag.1.percent+gdp.lag.1.1000+unemp.percent.lag.1+c_tf_r+dis_inc.lag.1, cluster = 'nuts_code2013')
t6_w5se<- coeftest(t6_within5, vcov=vcovHC(t6_within5, type="HC1")) 


#Table 6
t6<- stargazer::stargazer(t6_within1, t6_within2, t6_within3, t6_within4, t6_within5, 
se = list(t6_w1se[,2], t6_w2se[,2], t6_w3se[,2],t6_w4se[,2], t6_w5se[,2]), 
omit=c("nuts_code2013", "Constant", 'Year'), omit.stat = c('rsq', 'adj.rsq', 'ser', 'f'),
dep.var.caption = "Vote for Far-Right Parties", 
order = c("AGE", "BA_plus", "nojob", "emig.lag.1.percent", "imig.lag.1.percent", "gdp.lag.1.1000",
"unemp.percent.lag.1", "c_tf_r", "dis_inc.lag.1"),
covariate.labels = c("Age", "Education (BA)", "Unemployed","Emigration", "Immigration", "GDP", "Unemployment", "Current Transfers", "Disposable Income"), 
add.lines = list(c("Year FE", "\\checkmark", '\\checkmark',
'\\checkmark', '\\checkmark', '\\checkmark')))


